Three-dimensional theory of stimulated Raman scattering 
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We present a three-dimensional theory of stimulated Raman scattering (SRS) or superradiance. 
In particular we address how the spatial and temporal properties of the generated SRS beam, or 
Stokes beam, of radiation depends on the spatial properties of the gain medium. Maxwell equations 
for the Stokes field operators and of the atomic operators are solved analytically and a correlation 
function for the Stokes field is derived. In the analysis we identify a superradiating part of the 
Stokes radiation that exhibit beam characteristics. We show how the intensity in this beam builds 
up in time and at some point largely dominates the total Stokes radiation of the gain medium. We 
show how the SRS depends on geometric factors such as the Fresnel number and the optical depth, 
and that in fact these two factors are the only factors describing the coherent radiation. 



I. INTRODUCTION 

The collective emission of radiation from an ensemble 
of atoms is interesting both from fundamental as well as 
an applied perspective. The enhanced collective emission 
or superradiance [l[ from an atomic ensemble was pre- 
dicted already by Dicke in 1954 0] and was observed in 
the form stimulated Raman scattering (SRS) in 1962 [3|, 
but in recent years it has attracted renewed interest due 
to the the observation of SRS from Bose Einstein con- 
densates Ref. @, H, H, 0] • From a more applied perspec- 
tive the problem of SRS is closely related to free electron 
lasers Q as well as to activities in quantum information 
science aiming at realizing a quantum interface between 
light and matter Q . Recently it has even been proposed 
that SRS from a Bose-Einstein condensate could serve as 
a direct source of entanglement [l(| • 

From a theoretical perspective one of the challenges 
consists of describing SRS from an extended ensemble. 
Whereas the original Dicke superradiance [|J was de- 
scribed for a collection of atoms localized to dimensions 
much less than the wavelength of the outgoing light, 
most experiments are actually performed in the opposite 
regime where the dimensions of the ensemble is much 
larger than the wavelength. A full quantum descrip- 
tion of SRS was presented by Raymer and Mostowski 
in Ref. [ll[ using a one dimensional model. Such a one- 
dimensional description can be shown to be applicable to 
all transverse modes of the field if the gain medium has 
an infinite transverse extension [9J]. In such a descrip- 
tion there is, however, no restriction on the transverse 
modes and a summation over all transverse modes there- 
fore results in an infinite intensity of the outgoing light. 
The theory was generalized to also include certain three- 
dimensional properties of the propagation of light in the 
gain medium in Ref. [l2T |. Here it was argued that the 
one-dimensional theory could be used to predict the to- 
tal intensity for a sample with a Fresnel number of unity 
T = 1, where the process was dominated by a single 
transverse mode. These theories were developed under 
the basic assumption that the region in which this SRS 
process happens is defined by the properties of the laser 



both in time and space. Thus figures of merits are the 
width and temporal shape of the laser which is driving 
the SRS process. The experiments exploring the SRS 
process have changed since then @, H, d, 0|, and much 
more attention is given to systems where the temporal 
and spatial shape of the laser have long surpassed the 
spatial geometries and temporal properties of the gain 
medium. A three dimensional theory applicable to small 
atomic ensembles were presented in Refs. [3, [ljj in 
the approximation that certain off-diagonal matrix ele- 
ments in the momentum representation could be ignored. 
For the closely related problem of light emission from an 
ensemble of atoms with a few collective excitations, di- 
rect numerical simulations have been performed for a few 
thousand atoms [H, [l(| . Using periodic boundary con- 
ditions the three dimensional effects of this problem was 
studied in Ref. [17| , and an approximate analytical treat- 
ment was also presented in Ref. [18|. To our knowledge, 
however, no theory have been developed which fully de- 
scribe SRS from a spatially extended ensemble. Here 
we develop the theoretical framework that enables us to 
describe SRS from such extended ensembles. Since we 
shall neglect the depletion of the initial atomic state, the 
present theory is, however, only capable of describing the 
onset and build up of SRS. In the resulting theory the 
only two parameters describing superradiance is the op- 
tical depth d and the Fresnel number of the sample T '. 
We show explicitly that the time at which SRS begins to 
dominate is given almost exclusively by the optical depth 
but that the Fresnel number T is important for determin- 
ing the total amount of light radiated from the ensemble. 
Our theory is based on a generalization of the one dimen- 
sional theory presented in Ref. ll| , but includes several 



effects omitted in the three dimensional generalization in 
Ref. [HI]. The main difference is that we go beyond the 
extreme paraxial approximation used there, a generaliza- 
tion only briefly discussed in an appendix of Ref. [HI]. 
The theory that we develop can therefore explain both 
spontaneous emission as well as SRS. 

The analysis begins with the basic set of equations de- 
scribing the interaction of light with atoms. The atoms 
are treated as non-moving point particles and the radi- 
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ation fields are described by the displaced electric field, 
suited for a macroscopic description of the system. See 
e.g. Ref. for a discussion of this choice. We will 
in Sec. [Til derive effective equations of motion for both 
the radiation field and the atoms. These equations are 
directly comparable to the equations used in Ref. [ill ]. 
Having established the equations of motion we will in Sec. 
IIIII change from the point particle picture to a continuous 
description. This again follows methods described in e.g. 
Ref. [19l| . In Sec. II VI we make a formal diagonalization 
of the matrix describing the interaction between atoms 
mediated by the light. This diagonalization means that 
we have to find a basis that will simplify the interaction. 
In Sec. |V] we will look at the radiated field and see how 
this is evolving as the atoms are interacting. Finally in 
Sec. I VII we look at the intensity of the radiated field 
and present the final results. We shall in addition to 
the analytical results make a comparison with numerical 
calculations for the SRS starting with the point particle 
equations of motion derived in Sec. |TT] In Sec. I VIII we 
conclude the work. 
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Figure 1: Atomic level structure. Fwo stable ground states |1) 
and 1 2) are coupled through an exited state |3). We assume a 
strong classical laser of a + -polarized light drives the transition 
from |2) to |3) with detuning A. The laser thereby effectively 
drives a transition from level |2) to |1). The radiation uis 
connected to the transition from |3) to |1) describes the Stokes 
field, that is analyzed here. 



A. Atomic dynamics 



II. EQUATIONS OF MOTION 

In the electric dipole approximation the Hamiltonian 
describing a collection of atoms is given by 



H = J {H F +Hi}d 3 r + HA, 
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where D is the displaced electric field, B is the mag- 
netic field and P is the atomic polarization. The oper- 
ator a 3 nn — \n)(n\ is a projection operator for the j'th 
atom, and is the energy corresponding to the state 
\n). We choose to use the displaced electric field and 
not the electric field for reasons discussed e.g. in Ref. 
[lj|. This choice, however, does not influence the result 
of the analysis. Here we have ignored any direct inter- 
action between the atoms, e.g. atomic collisions. As we 
shall often make reference to Ref. [ll|, we shall try and 
match the constants and the dynamics of our system to 
the system presented there. The Hamiltonian is also cho- 
sen such that results derived in Ref. 1!) can be directly 
incorporated. In the following section we will focus on 
the dynamics of the atoms. 



The macroscopic description of the atomic ensemble is 
given by the polarization, P(r, t) which again is the sum 
of the individual dipole moment of the atoms. 



P(r,t) 



S(r - rj)d nm a J nm (t), 



(11.5) 



where the time dependent operator <J 3 nm (t) is the opera- 
tor |n)(m| taking the j'th atom from state |m) to state 
|n), and the dipole moment is d nrn = e(n|r|m). In ad- 
dition we assume the atoms to be identical with a level 
structure shown in Fig. [TJ We assume the two levels 
|1) and |2) to be stable ground states. For the chosen 
atomic system we assume that the transition from level 
|1) or |2) to |3) increases the atomic angular momentum 
by one unit of h, and that there are no other states that 
the level |3) can decay to. This means that the only non- 
vanishing vector components of the dipole moments are 
e + = (e^ + ie l) )/v / 2 for positively oscillating terms and 
el for negatively oscillating terms. 

We employ the Rotating Wave Approximation (RWA) 
and assume that A is sufficiently large so that we may 
adiabatically eliminate the exited level |3). In this pro- 
cess we split the radiation field D into its positively and 
negatively oscillating parts, and extract the strong clas- 
sical field V c i oscillating with a frequency ui L from the 
weak quantum mechanical stokes field D oscillating with 
frequency lu s . We will assume that the strong classical 
field is constant over the region of the atoms and can 
be written as a plane wave with a constant amplitude 



V 



(+) 



e + . The presence of the strong 
classical field V c i induce a Stark shift of the atomic lev- 
els. The effective Stokes frequency lu s is therefore given 
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by 



L0 S = W L + U) 2 1 



\d 3 i\ 2 \Vd 
7i 2 e 2 A 



(11.6) 



We define slowly oscillating operators both for the atomic 
operator 021 and for the stokes field D 



(11.7) 
(II.8) 



For large detuning and weak fields we can adiabatically 
eliminate the exited state, and obtain an effective ground 
state equation of motion. 

Jt^ {t) = 5 (a ' 22 _ ^WDtir^t), (II.9) 



where the constant a is given by 
Cl ~ he A ' 



(11.10) 



The positively oscillating part of the polarization is in 
this approximation 

P (+) (M) - E ^d\B + a[ 2 (t)5(v rj). (11.11) 



The negatively oscillating part p( \r, t) is found by Her- 
mitian conjugation. 



B. Field equation 

The equation of motion for the electric field D(r, t) is 
e.g. given in Ref. [l9j by 



D«(r,t)=D+(r,i)- 



dt' pW(r, t\rj, t') ■ e+a\V cl \d( 2 (t') 



(11.12) 



where Do is the unperturbed field containing the vac- 
uum Stokes field and the classical laser-field, and P^ 
is the propagator. The coupling between level |2) and 
1 3) in principle giv e rise to an index of refraction. As 
shown in Ref. [191 ]. such an index of refraction should 
be incorporated into the propagator P^K In the limit 
of large detuning A (but fixed a|2? c /|), we can however 
neglect this, and will do so in the following. The prop- 
agator in the slowly varying approximation is in Fourier 
representation given by 

r . . 1.2 ik'(r-r') 

p(+)(r,r') = I rf3 fc ^ _^£_ ix ££% (n.13) 



(2tt) 3 (/c 2 - 1) 



where the k-integral is understood to include only the 
contribution corresponding to the retarded Green func- 
tion. Here and in the remainder of this work we will mea- 
sure the spatial coordinates in units of k s , which gives the 
factor of and a pole at 1 in Eq. (|II.13[) . 

Inserting Eq. (III.12j) into Eq. ()II.9[) gives us an effec- 
tive equation of motion for the atomic operators, 

jfiS) = -~4>(*) + E Mii'tkt) + Fj '(*), (H.14) 



where 
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F 3 (t) 



Za DW(r,,i).e;^'(t), 
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e* + -P(+\r j ,r' j )-e + . 



(11.15) 
(11.16) 
(11.17) 



We have in addition made the approximation 022 — 011 ~ 
1, where we assume that initially all atoms are in state |2) 
and that the experiment takes place on a timescale such 
that we may neglect depletion of this level. To derive the 
decay T we used the identity 



el.p(+)(r,r J ) 



67T ' 



(11.18) 



which is discussed in e.g. Ref. [19( as the infinitely short 
propagator, and the relation (022 — &\\)<J\i = 
The effective equation of motion for the atoms. (III. 14p 
is the starting point for many studies of SRS [ill Il5l|. 
but also for studies of the coupling between atomic spin- 
excitations and collective emission of light, [HI [lj], l20| . 
In our analysis we neglect the effect of the source term 
Fj in Eq. (|II.14[) , as we are eventually only interested in 
measuring the ph oton flux !>(+))■ It can be found 

from Eqs. pi,12j) and (|II.14j) that the effect of the source 

term Fj leads to a contribution (Dq ' Dq + ^) to the mea- 
surement. This contribution vanish as we assume that 
the Stokes field is in the vacuum state. We also assume 
that there is no classical noise in the laser field 2? c ;. 

We shall be interested in defining creation and annihi- 
lation operators for the atoms. This leads in general to 
nonlinear equations, but under the low excitation approx- 



imation. 



that is (T22 



a 1± ss 1, we employ the Holstein- 



Primakoff approximation and simply use 



b) =0-12, 



so that 



bj = (T21 



(11.19) 



(11.20) 



The effective equation of motion for the atoms is then 
given by 



d = + E Mj/b],®, (11.21) 



dt 
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and for the field Eq. (|II.12|) gives 

D<+>(r,t)=D+(r,i)+ 

I dt' P^\v,t\r^t') ■ e + a\V cl \b]{t). (11.22) 



III. GOING FROM DISCRETE TO 
CONTINUOUS SYSTEM 

We will be interested in treating Eq. (|II.21|) as a con- 
tinuous equation. For an atomic gas we do not know the 
individual positions of the atoms, thus an expectation 
value of a physical operator has to be accompanied by 
a spatial average of the individual atomic positions. We 
therefore define the density distribution /5(r), 

p(r)=^5(r-r,-). (HU) 

3 

We assume that after a spatial average of the position of 
the atoms in the ensemble the density distribution p(r) 
can be described by a Gaussian function 

(p(r)) sa . =p(r)= Po e ^ ^u. (III.2) 

We will also assume that 1 <C a± <C <j\\ and a\ > &\\ 
where spatial coordinates are measured in units of k s . 
We then define the normalized continuous operator 

b(T) = -jL=Y,8(T-T j )b j . (IIL3) 

After taking spatial average of the position of the atoms, 
this definition leads to the standard commutation rela- 
tions for such continuous operators, 

[b(r),P(r)] = 5{r-r'). (III.4) 

From this definition of the continuous operators Eq. 
(|II.21[) can be rewritten 

iUt(r,t) = / d\ Y^ 5 -^^M{v,v')^)b^v\t) 

= J d\ ^)M(v,r')y/p~(?)b\r\t) 

+ f rfV S{T - r ^Z P{r) M(r,r')^(?)bHr',t). 
J j VP( r ) 

(111.5) 

The lowest order spatial average is found simply by mak- 
ing a spatial average of Eq. (|III.5|) . In Ref. [19( we 
considered higher order corrections coming from such a 
spatial average, and showed how fluctuations in position 
give rise to spontaneous emission and dipole-dipole in- 
teraction effects. Here we shall ignore these effects. To 



lowest order in the spatial average, the first term in Eq. 
pil.5l) describes the mean effect of the atoms interac- 
tion with each other, that is when averaged with respect 
to their individual positions. The second term will af- 
ter spatial averaging only give a contribution for atoms 
interacting with themselves via the infinitely short prop- 
agator 19], thus the term effects in the decay described 
by r, which is independent of the interactions between 
atoms. To get the sign of the decay, one would have to 
remember that the approximation er 2 2 — en f» 1 is not 
justified for this particular type of term, and including 
this correction as in Eq. (|II.14[) . gives the negative sign. 
The continuous version of Eq. (|II.21| is then 

j/{r,t) = J d?r ^plv)M{v,v>)^pl?)b\v> ,t) 

-\b\v,t). (III.6) 

It is convenient to remove the last term of Eq. (|III.6|) by 
defining new atomic operators with respect to the decay 
T, [ w(r,t) — ► ¥(r, i)e~^ 1 ], ignoring the source term 
F and the point-particle corrections, the effective differ- 
ential equation describing the excitation of the atoms is 
after spatial average given by 

j/(r,t) = J d 3 r V^Af(r,r')V^>V ',*). 

(III.7) 

Similarly the field equation (|II.12[) can be described in 
terms of the continuous operators, and one find 

DW(r,i)=D+(r,t)+ 

a\V cl \ J d 3 r'pW(r,r') v / Pp)-e+6 t (r',i). 

(HI.8) 

In the following we will find approximate solutions to the 
above equations. 

IV. DIAGONALIZING THE INTERACTION 
MATRIX 

The system is assumed to be cylindrically symmetric, 
with a density described by Eq. (|III.2[) . We shall there- 
fore use a cylindrically symmetric set of basis functions 
for our diagonalization: a combination of plane waves 
and Bessel functions. We denote the basis by {fkmn}, 
where 

\/2 r 

f, (r Z (h\ = - ikz+irrup r l v \ 

ZTra c J m+ i{A. mn ) a c 

(IV.l) 

J m is the Bessel function of first kind of order m, and 
X mn is the n'th zero of the m'th order Bessel function 
of first kind. The parameter a c is a cut-off in the radial 
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Figure 2: Sketch of the integration contour Ci, in the integral 
representation ()IV.5|I of the Green function. 



direction, meaning that our basis is complete on the in- 
terval r € [0, a c ]. The inner product defined for this basis 
is therefore given by, 



dz / rdr 9*{r 1 z, (j))ip(r, z, c 
oo Jo 



(IV.2) 



For a discussion of this basis see e.g. Ref. [21(. To solve 
Eq. pil.7|) we will diagonalize the matrix given by 



M(r,r') = -^^e] 



v ^)p(+)(r,r')VpM.e + . 

(IV.3) 



The propagator P(+) is found in a real space represen- 
tation in e.g. Ref. [22]. One may from the real space 
representation of the propagator show that 



e; • P(+) (r, r') • e + = _*£ (V 2 + df) t^l 



8n 



r — r 



(IV.4) 



The polarization effects are here included in the differen- 
tial operator V 2 + (9 2 . In addition we use that the Green's 
function can be written as [291 



i\r— r | 



_ — iST f r ] hp im(4>-4>')+th(z-z') 



J m ( Vl - h? T< )H£){ y/l - h 2 r> ) , 

(IV.5) 



where r< (r>) is the minor (larger) of r and r'. C± is 
describing a curve essentially going from — oo to oo along 
the real axis but shifted to avoid the branch cut and pick 
out the retarded Green's function, as shown in Fig. O By 
introducing an integral, the non-trivial product of Bessel 
functions in Eq. (|FV.5j) . can be symmetrized |24l |: 



J m ( VI - h 2 r < )H W ( - h 2 r> ) 



in J x 2 + h 2 — 1 



(IV.6) 



The propagator is then given by 



[ .'"'I 1 + 1,2 



/l 2 - 1 



e im( -' t '-^ +ih( - z - z ' ) J m (xr)J m {xr'). 

(IV. 7) 



In the basis {fkmn} the differential equation (|III.7|) can 
be written 



dt 



ft (f\ — \^ f.^kmn it 
u kmn\ L ) ~ /-^ lvl k'm'n' u k'm'n' 



(IV:, 



k' m'n' 



where 



=</ fcmn (r)|M(r,r')|/ fcw (r')), (IV.9) 



and 



&Ln(*)=</*mn(r)|6t(r,t)). 



(IV.10) 



When calculating the matrix Eq. pV.9[) . we have to make 
integrals over r, z and <fr. We can at this point simplify 
the radial integrals by extending the upper integral limit 
to infinity. This is correct since the cut-off a c can be cho- 
sen arbitrarily and as we in the end will set it to infinity. 
Due to finite width a± of the density function, this limit 
accurately describe the matrix elements for a c 3> <t_l ■ Af- 
ter making the spatial integrations the matrix M reduces 
to 



M£"=Wy dh xdxr,(k-h) V (k'-h)- 



I m (2a\ ln x)I m {^o\ ln ,x) (IV. 11) 
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where 

r)(k) =%e- a l k \ (IV.12) 

and where we have introduced the constant Ao = 37f o° ■ 
To shorten notation we also introduce 7„ = TO " , where 
we understand that j n depend on the azimuthal quantum 
number to. For integrals involving Gaussian functions 
and Bessel functions we refer to Ref. [24j. We notice 
that both integrals over x and h are bounded by Gaussian 
functions, and since we assume a± 3> 1 we may make a 



series expansion in x and h of the function l/(x 2 + h 2 — 
1). We will be interested in a series expansion of the 
integrals over x and h only to the lowest order. Since 
we assume that an ^ erj_, i.e. cigar-shape, our lowest 
order calculation will terminate after first order in 1 ja\ . 
The integral over h can to this order be approximated by 
treating the function rj(k — h) as a delta function, thus 
we shall here and in the remainder of the article treat 
the function r/(k — h') as a delta function. We show in 
Appendix [X] that the integral over x to lowest order in 
the variable l/oj_ gives 



M-, 



k m'W 



An 



A r 



A 



1 1 n 



k 2 - 1 VSa 2 ± {k 2 - 1) 



O 



(IV.13) 



where 



A" 



(IV. 14) 



and 



4 2 7„+7„/ r 2 ( t 2 'y 7 / 
Al m _ 4g J- e i "4 z "J-7«7n'J mr ic-i 

" c Jfn-\-l y^-mn \^mn' ) 

The matrices A"i, and A 1 ™ / are normalized such that 

fill I LI L 

for a± — > oo they reduce to a delta- function <5(n — n'). 

In the following we take a closer look at the matrix 
A" l n , defined in Eq. (|IV.14|) . For simplicity we will not 
consider the correction A 1 nn /, however the conclusions 
drawn in the following holds for the correction as well. 
The differential equation for our system with respect to 
the quantum number n, n! has got the form 

j t b n (t) = J2^Kn'K'(t), (IV.16) 

n' 

where u> is some real number. We wish to take the limit 
a c — * oo. To clarify what this means let us write the 
matrix A in the following way: 

_2„2 



Am _ a i, ~m 2 2 -- 



(IV. 17) 



where 



1—1 nri/ 



i 



^ V ' -X-mnXnin' Jm~\-1 (^mn ) Jm-\-\ \P^mn' ) 

(IV. 18a) 

(-!)"+"' for A mn ,A„ m ^oo 

^ (IV. 18b) 



^■kjnn' — 

a c 



(IV. 18c) 



We thus see that when letting a c — ► oo, a transverse 
momentum naturally arises fc m j_ = lini ac ^oo k mn , and 
the discrete matrix equation, Eq. (|IV16[) becomes an 
integral equation over the transverse momentum fc m j_, 
using J2 n > Afc m „/ -> / dk m± 

^b(k m± t) = J dk' m± nA m (k m± ,k' m± )b(k' m± ,t)- 

(IV.19) 

It is evident that when using the limiting properties of the 
Bessel function I m (x) the integral kernel A m (k m ± 7 k' ml ) 
becomes a delta function for a± — > oo. 

7T 2 <ri e - ^ +fc ' ™x ) J m (tt 2 ^ fcm± fc ^ ) ^ 



(IV.20) 

We thus realize that the effective one-dimensional result 
obtained by Raymer and Mostowski [ll[ is exact for ev- 
ery transverse mode in an infinitely wide atomic ensem- 
ble. In this limit however there is no limitations on the 
transverse momentum, which results in an infinite inten- 
sity. To obtain finite results we thus need to consider the 
full solution to the three dimensional problem. 

Now we again include the correction A 1 nn , in the anal- 
ysis. Both matrices A™ n , and A 1 ™^, are real and sym- 
metric and can thus be diagonalized. In Appendix [B] 
we show that the two matrices commute. We can there- 
fore choose a common set of eigenfunctions, |-F)tinn( r *)| 

for both matrices. We define the unitary matrix U that 
transform our initial basis {fkmp} to the basis given by 

the eigenfunctions |i 7 fe m n(?')}i 



(IV.21) 
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Finally we will define a corresponding set of eigenvalues, 
A pp' = £ U\ n \ mn U np , (IV.22) 



and 



A pp' ~~ ^ ] p n \ mn U, 



np' 



(IV.23) 



It is convenient in the following to change to this basis, 
where A™ n , and A v ™ n , are diagonal. We therefore write 
Eq. (1TVT13T> as 



Errt M k ' m ' p 'TI , , - 

PP' 

1 4- k A mn 



1 + fc 2 



o 



„-2 „-4 
0|| 7 J_ 



(IV.24) 



where the index s refers to being at the symmetry axis. 
To arrive at the above result we used the real space repre- 
sentation of the propagator P Eq. to leading order 
in one over distance. This approximation is done out of 
convenience but is not strictly necessary. When calculat- 
ing the mode expansion of the electric field in the general 
modes Fkmn we shall check that the limit t — > exist and 
is given by the expression, (jV.lj) . 

The analysis of the radiation field for t ^ starts by 
inserting the identity operator, 



1 = / d 3 r' J dk^F kmn (r)FZ mn (r>) 



(V.2) 



into the field equation, (|III.8p . We then get the following 
expansion of the electric field. 



V. REAL SPACE REPRESENTATION OF THE 
ELECTRIC FIELD 

In the following section we will, based on the eigenvalue 
analysis of the atomic operators, derive the real-space be- 
havior of the electric field. We shall divide the analysis 
into a regime of small times where the dominating effect 
is spontaneous emission, and a large time regime, where 
the dominating effect is the cooperatively emitted light, 
the SRS beam. To keep things simple, we mainly con- 
sider the electric field at and around the symmetry axis. 
In this region the scattered radiation field is sufficiently 

well described by the vector component and its Her- 
mitian conjugate. This can be seen from Eq. (]III.8[) and 
the real space representation of the propagator (|II.13|> . 

Let us first determine the electric field on the symmetry 
axis at the initial time, t — 0. In this case the electric 
field is given by: 



^ +) (r s ,0) ^ +) (r s ,0) 



d 3 r 



, a\V cl \kl 



4tt 



&t(r',0): 



(r 2 + {z- Z ') 2 f/ 2 



(V.l) 



D^ + (r,t)=D ( +\r,t) 

+ f dV f dk^kmn^e^F^'p (r',0), 

(V.3) 



where 



C kmn {v) =a\V cl \ j d 3 r'e;-#( + )(r,r')-e +v ^OF fcm „(r'), 

(V.4) 



the functions Fkmn are the basis functions given in Eq. 
(|IV.21j) , and the eigenvalue Xkmn is given in Eq. (|IV.24|) . 

The calculation of the modefunctions Ckmn is initiated 
by integrating with respect to the spatial coordinate r'. 
The integrals involving Bessel functions are found in e.g. 
Ref. H, and one arrive at 



Cftmn(r) = - / dy / xdx }_^U np — — — -x 

47r J J x l + (k + yY - 1 a c J m+ i(X mp ) 



I — 



e -° n v -'xl7,+x n m {2ai lp x). (V.5) 



The next step of the calculation is to include the mode summation. We will therefore define the propagator p(+) 



8 



given by 



pW(r,r';t) = J 



dk 



(V.6) 

We notice that the variable y in Eq. (|V.5|) is small, as 
it is controlled by the Gaussian function of width 1/cti i . 
We shall therefore by a translation of the integral vari- 
able k' = k + y move the perturbation y to the eigen- 
value Xkmm so that we use Xk'- y ,mn- This choice ensure 
that we will get the correct behavior of the integrals in 
the limit t = 0. By doing this we can then in principle 
make the k' integral by using the series expansion of the 
function e Afc '-»< m ™*, where the zeroth order term in the 
expansion in t is the limit given by Eq. (|V.1[) . In order to 
accurately capture the exponential growth, we however, 
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instead follow the path used by e.g. Ref. 

In the following we make a series expansion of the 
eigenvalue \k-y,mn given in Eq. (|IV.24[) with respect 
to the variable y. 



A 



k—y,mn 



= -T I K. 



{k-yf + l \l n (k-yf + 1 



l {k-yf-l V8 CT 2 ((fc _ y) 2 _ 1)2 
fc 2 + l 



1 



'x 
k 2 + l 



(A- 2 



1) 



(V.7) 



The series expansion can be done since the y-integral is 
bounded by a Gaussian function. To shorten notation we 
have substituted k' — > k, and introduced the coefficient 
- \ _ „ 

In Eq. (IV. 6p the fc-integral includes a pole 
1 1 



k 2 + x 2 - 1 



(V.8) 



2VT~^ a?(fc - v 7 ! 3 ^ 2 ) 

where the arrow reflects the fact that we are only in- 
terested in the retarded Green function, which corre- 
spond to the pole k = x 2 . Since we are particu- 
larly interested in this pole, we shall in the fc-integral 
in Eq. (|V.6|) . make a translation of the eigenvalue 
1 _ x i mn , and then a series expansion 



k—y,mn 



A 



k-y+V 



similar to Eq. (|V.7[) . We can make the calculation with 
two different situations in mind: One situation explains 
the spontaneous radiation originating from a sample of 
atoms of some geometrical shape. We are most interested 
in the other situation describing the collective emission 
or the SRS occurring when the atoms co-radiate. As a 
check of our formalism we shall, however, also consider 
the short time-limit where there is just spontaneous emis- 
sion. We expect that as time evolves the SRS effect will 
become dominant. Therefore we demonstrate where the 
SRS effect is found and described in our mathematical 
treatment of the problem. 

Let us first show how the important steps in the calcu- 
lation of SRS is done, before going into the full details. 
The integral appearing in the calculation is of the type 

o^k-v ,mnt+ikAz 



21t(t) 



1 

2^ 



dk- 



k 2 



1 



(V.9) 



where Az — z — z' . For now we consider the lowest or- 
der correction for simplicity, that is we neglect /i m „ in 
Eq. (|V.7jl . Including \i mn to the eigenvalue is a trivial 
generalization. We focus on the pole in the integral at 
k = Vl — x 2 , as this pole describes the energetically al- 
lowed scattering processes. By introducing the variable 
s = iAz(k — yl — x 2 ) the integral l£ can be written 



1 



ds- 



2V1 



(V.10) 

where the superscript indicates that this is a zeroth or- 
der calculation in the correction to the eigenvalue due 
to finite size. The SRS contribution to Eq. (IV.lOj) 
comes from the pole of the exponential. In order for 
this pole to contribute to the pole describing the prop- 
agated light, that is the zero point of the denomina- 



tor, the term Az(vl 



Az(VT 



1) has to be small. For 



1) < 1 we shall treat it as a perturbation. 



When this no longer apply, the pole in the exponent can 
be neglected, and we are thus left with the result for short 
times, i.e. spontaneous emission. The latter is analyzed 
in the following section, and we shall for now concern our- 
selves with the SRS contribution. For reasons discussed 
in Sec. IVBI wc will, when discussing SRS, use that Az is 



large, so that Az(Vl — x 2 — 1) 



2 A 2 



Since 



x Az 



2 • i 2 

we can make an expansion in this quantity and obtain 



< 1 



T k {t) 



1 



-2itfi 



ds- 



s l+l+2q 



(V.ll) 



Here we include the correction to the eigenvale in Eq. 
(|V.7[) . The integral may be found in Ref. [25[ and we 
find 



A A ? OO OO 



= VEE 



1=0 q=0 



2 Az\ l (^awAz 2 ) 9 



Il+2q{2V KnntAz) 

y\ mn tAzy+^- 



(V.12) 



A. Short time limit 

In order to understand our calculation of SRS, we first 
analyze it for t = 0, as we know how the propagator for 
t = looks when measured on the symmetry axis. The 

2 A 

t = regime is also met for x ^ > 1 . We shall also refer 
to this calculation as the short time limit. Here we find 
from a residue calculation Eq. (|V.9|) to give 

ipWl — x' 2 Az 

r ' (0) " -^T=W (V13) 



By inserting this into the propagator in Eq. (|V.6|) . the 
propagator may be written 
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P«(r,r';0)=£ 



a\V cl \kl^p- {l-dl 

47T 



np' ' 



J m {xr)J m { lp ,r')I m {2a\ lp x) -*Ai p +* )-^r 

i 



(V.14) 



where A0 = <f> — <f>' . The only dependence on the mode- 
index n is in the product of the two matrices U np U np i 
and the sum over n reduces to a delta function 5 PP ' . We 

then, similar to Sec. IIV1 identify ^- — * J for 
a c — > oo. The variable j n is in this sense fixed, thus 
letting a c — > oo has to be accompanied by -X" mn — > oo. 
Therefore we can use the large argument approximation 
for the Bessel functions, 



Using this we can make the integrals over j p and j p i 
The result of the mode summation (|V.6|) is then 



2 m?r 7T 

— — cos(X mn -_--), X 

7TA mn A 4 



(V.15) 



J 



P ( + ) (r y. 0) _ |P c ,|fc3Vp(P)(l-^) Tc „ nl 



<s- 



^gi\/l-£c 2 Az 



" n , J m (xr)J m (xr') (V.16) 

VI — a; 



This is the main result of this section. To verify the 
validity of the approach taken so far, we shall now show 
that the propagator (|V.16|) reduces to the one found on 
the symmetry axis, (IV.1[) . In order to show this we will 
use the summation theorem for Bessel functions, see e.g. 

a 

Y / e mA *J m (xr)J m (xr r ) = J (xR), (V.17) 



where R = W r 2 + r' 2 — 2rr' cos(A0). In this way the 
propagator in Eq. (|V.16|) can be written 



P ( + ) (r y )0) _ a\v cl \kl^{7)(i~dl) y _ 



^i^/l — x 2 Az 



cdx- 



MxR) 



(V.18) 



The ^-integral is known and may be found in Ref. [2o| . 
to give 



P (+) (r,r',0) = 



-a\V cl \k 3 s ^(F)(l - d 2 ) e ^+^ 



8tt 



VP 2 + Az 2 
(V.19) 



Finally the z differential give us the result we are looking 
for. 



P (+) (r,r',0) = 



a\V c i\kl^(?) e'v^+Ai 7 ±R2 + Az 2 

4^ VP 2 + Az 2 R 2 + Az 2 ' 

(V.20) 



When we then look at the symmetry axis, the variable R 
reduce to r' and we are left with the result in Eq. (|V.1|) . 
The result of this section can be written as 

D { +\r,0) =L^ +) (r,0) o + J dV P« (r, r'; 0)^ (r', 0). 

(V.21) 

B. Finite time, build up of SRS 

In the following we shall analyze the effect of the eigen- 
values X mn and Xl nn in the expression (fVTT2|) . When we 
introduced the eigenvalues in Sec. IIVI we only concluded 
they could be found. We also know that physics con- 
nected to the eigenvalues can not depend on the cut-off 
a c involved in the index n. In the following we show that 
indeed the physics is independent of the cut-off a c . To 
find this result we shall in particular look at the sum 
Z~2 n ^ pnVmn^mnUnp' where the powers N and M are 
zero or some positive integer. [ The powers N and M are 
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connected to the series expansions of functions involving 
the eigenvalue A m „, e.g. Eq. (|V.12|) . ] \ mn and X$ nn are 
the eigenvalues of the matrices A™, and h} pp , in Eqs. 
(IIV.22P and pV.23p . Let us generalize the matrices A™, 
and A 1 ™, defined in Eqs. PV.14P and PV.15P to 



vf\ N 



4a ± 2 e ^^ 2+ ^ 2 )/ m (^ 7p7 p') 



TV J m(-^mp) Jm {^mp 1 ) 



(V.22) 



i.e. A™, correspond to N = 2 and A pp/ correspond to 



iV = 1. One can then show that 



This result along with the appropriate series expansion 
of functions involving the eigenvalues X mn and \ mn can 
be inserted into the result for the propagator Eq. (|V.6|1 . 
and the resulting sum over indices p and p' takes the form 



- J m {jp>r')l m (2a± 2 7px)e- 



pp' 



' Jm+1 (A mp )t/ m _(_x 



4<7_L 



(X mp i) 
Jm {XT ) j 



pp'\ N ) 



(V.24) 



M , v 



2(JV + M — s) + s 



(V.23) 



where N is an integer derived from Eq. (|V.23p and 
the before mentioned series expansions. The propagator 
(|V.6P can therefore be written 



p (+) (rj r ':t)=J2 a ^ff~ d ^ I xdx I dy^Le-^^' £ U np U np , x 

mn * pp' 

, 2 imA^ T ■/m(gT-)Jm(7 P '^)^m(2(7i7pig) - a *^+ x ')-£ 

4o±e ^fc— J7 7^ — w 7^ T 

t*c ^mp J^m+l l^-mp' J 

w|2? ci |fc s 3 



4?r 



V^)^ / ' xdxe mA ^ +lAz J m (xr)J m (xr')x 



oo oo 



EE ^ ^ r '-''» /U v *» (v ' 25) 



where 



We notice since (x 2 Az)/2 < 1, that choosing the vari- 
able Az large means that the sum over I will converge 
very fast. Choosing the variable Az large can be done by 
placing the detector plane far away from the sample, in 
which case we will talk about a far-field calculation. Un- 
fortunately the sum over q converges more slowly when 
Az is larger, and we can not quite rely on our initial ap- 
proximations [r)(k — k') w 5(k — k'), see Sec. |IV] for large 
Az. We shall therefore consider the problem in the near 
field region. The limit ^2/Az in the x-integral we shall 
on the other hand approximate with the value y/2/L, 
where L = \/2waii is the effective length of the atomic 



ensemble. This approximation will become better at later 
times, since the coherent build-up is essentially described 
by the modified Bessel function Ii+2q{2\J XaAzt) which in 
time will dominate for large values of Az. In Fig. [3] we 
illustrate the physical significance of the integral over x, 
which represents an integral over transverse momentum. 
We see that as we include more light from deviating an- 
gles, this radiation has a shorter region over which it 
can build up, and as the build-up is exponential in the 
build-up length, the error made by the cut-off L becomes 
relatively small. From the propagator (|V.25P the electric 
field can be written, similar to the spontaneously emitted 
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Figure 3: A sketch of the coherent build-up of radiation in 
an atomic cloud. In principle the build-up can happen along 
any direction, however for a cigar-shaped geometry the most 
significant build-up happens along the axis of the cigar. 



radiation, (|V.21|) . as 

D [ +\r,t) = D ( | +) (r,i) + / dV P<+) (r, r'; 0{r', 0). 

(V.27) 

VI. INTENSITY AND THE CORRELATION 
FUNCTION 

In this section we consider the electric field, and as- 
sume that we place a detector in a plane at some position 
zq after the end of the atomic sample. We then define the 
correlation function as a function of the radial coordinate 
r and time t 

C ( r > r '< *) = f <MD+H*o, r, 0, t)b { + ] (z , r', <j>, t)>, 
neoh J 

(VI.l) 



where (•) is the quantum mechanical average. The nor- 
malization g 2 , is chosen such that the number of pho- 
tons in a pulse is given by 

N p = JjgJ dtC(r,r,t). (VI.2) 

The factor fef is inserted since lengths are measured in 
units of fc s . Inserting the propagator in Eq. (|V.25[) allows 
us to describe SRS, while the propagator (|V.14[) gives the 
spontaneous emission for short times. We shall be most 
interested in SRS, but will also for comparison examine 
the spontaneously emitted light. First we present the 
correlation function describing the SRS, when measured 
in a plane at the end of the atomic sample. An important 
parameter below will be the Fresnel number T which we 

2 

define by T = -jK [Recall that all lengths are measured 
in units of k s .] We shall in general assume the Fresnel 
number to be large, in particular T > 1. In the integra- 
tion over z' we will use the following substitution 



J dz'e^n -> dz', (VI.3) 

where L = \f*hta\ \ . The correlation function can then be 
calculated to give 



where 



m Iqk n,n r 
I q k 



Xlqkn k\k'\{l + 2q + k)\{V + 2q> + k')\ \ ' 



E(n/2), n -2s, _ ( (-1) Q + Q ' + S+S ' (2tt)~ 



E(n'/2) n '-2s' 2 I -n')\(n-2s-Q)\(n' -2s' -Q')\s\s'\Q\Q'\ , 

lqkn ^—f {l+Q+Q' + k+k'+l+l'+2(q+q'))(2+2{k+k')+q+q'+n+n')' 



(VI.4) 



This is the main result of this section. We notice that 
when r is measured in units of a±, the only variables 
controlling the behavior of the correlation function is the 
Fresnel number, T , the optical depth, d — 6i:poL and 
time measured in units of the single atom scattering rate 



r. This follows since X ptL = ^-p LTt = From 
the correlation function (|VI.4[) we also expect fast con- 
vergence in the index q and I as the Fresnel number in- 
creases. In the remainder of this article we shall evaluate 
the correlation function numerically. Even though the 



12 



correlation function involves a double integral beside the 
large number of sums, we see that as we increase the 
index fc, g, q' , n, n', the y- and y'-integrals will sim- 
plify. This follows since the argument of the modified 



Bessel function decreases as the indices k, k', q, q' , n, n in- 
creases. We can therefore use the small argument limit. 
Similarly the Gaussian function can be approximated by 
unity. From Eq. (|VI.4|) we see that the dominating term 
in the sum over k will have a higher k when time grows. 
This means that the radial behavior of the beam simpli- 
fies. Due to the small argument description of the modi- 
fied Bessel function the radiation is eventually dominated 
by the m = mode. 




A. Intensity on the symmetry axis 

In this section we will examine the radiated light on the 
symmetry axis. The purpose is to examine the timescale 
on which there is a crossover from spontaneous emission 
to SRS. 

Placing the detector on the symmetry axis is a nice 
simplification especially for the spontaneous emission 
correlation function, since in that case we may use the re- 
sult presented in Eqs. (fV\2T|) and (|V?20| . Also the SRS 
correlation function simplifies since terms with m =/= 
vanish at the symmetry axis. In the spontaneous emis- 
sion limit t ss the intensity on the axis is given by 



Co(0,0) = fc s 2 A 



dAz / r'dr' . 



Az 2 ) 5 



(r' 2 



Az 2 ) 3 ' 
(VI.6) 



where we use the substitution in Eq. (|VI.3[) . and assume 
that the detector is placed at the end of the atomic en- 
semble. The z-integral can be performed analytically and 
one finds 



Figure 4: Plot of the time r c measured in units of (dF) _1 , 
at which the intensity on the symmetry axis is dominated by 
SRS. The cross-over time is only weakly dependent on Fresnel 
number, and is given primarily by the optical depth. 



optical depth. Thus in order to compare the two time 
domains, the spontaneous emission and the SRS, we will 
have to fix e.g. the length of the system. From Eq. 
(|VI.4p we find that the cross-over time when going from 
spontaneous emission to SRS scales linearly with the op- 
tical depth, so that an increase of the optical depth gives 
a similar decrease of the cross-over time. In Fig. |4] we 
show this cross-over for varying Fresnel numbers T and 
a fixed length of the ensemble L = We see that 

the cross-over only depends weakly on the Fresnel num- 
ber. The main parameter characterizing the time scale is 
thus the optical depth. The cross-over time is found by 
plotting the intensity on the symmetry axis, Eq. (|VI.4|) 
and the spontaneous emission on the symmetry axis, Eq. 
(|VI.7[) . and finding the point at which they cross. 



C Q (0,0) =fc 2 A L / rdr 



32 



-13-llr 2 19arctan(r _1 ) 



(l + r 3)2 



(VI- 7) 



From this expression we find that the parameters con- 
trolling the intensity on the symmetry axis is the optical 
depth, and the ratio between the length and the width 
of the atomic ensemble. 

We shall now investigate the time scale on which SRS 
begins to dominate the radiation. For short times where 
the radiation is dominated by spontaneous emission, we 
expect that the radiation is being emitted almost ho- 
mogeneously in all directions, so that it makes sence to 
compare the spontaneous emission in a given direction, 
with SRS. We find from Eq. (jVLT)) that the figure of 
merit for the spontaneous emission is the density, the 
length, and the width of the atomic ensemble, and not 
as in the case of SRS, only the Fresnel number and the 



B. Intensity profile 

In this section we shall look at the spatial shape of the 
radiation leaving the atomic ensemble. Before we present 
the numerical calculations for the coherent emission we 
will look at the correlation function in Eq. (|VI.4[) . The 
spatial shape of the function is mainly given by 



2jr f 2:F r ,-r' 

dy / dy'J m (^/y — )J m {\nj l — )x 
o Jo f-l CT -L 



e 2 + 2(fc + fc') + 9 + <!' + n + n') J 



2+2(k+fc') + g+g' + " + ™' 



(VI.8) 



With increasing values of k,k' , q,q' , n and n' , the expo- 
nential function can to a higher and higher precision be 
approximated by unity. The modified Bessel function of 
order m can for small arguments be approximated with 
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an m'th order polynomial 



Im(z) 



(z/2)> 



ml 



(VI.9) 



From the argument of the modified Bessel function in Eq. 
(IVI.8I) we find that the region for which the approxima- 
tion Eq. (|VI.9|) is applicable is given both by the number 
2 + 2(k + k') + q + q' + n + n' and by the integration 
range IT . Eq. (|VI.8[) indicates that as time increases 
the dominant mode will be the to = mode for a finite 
sized atomic ensemble. On the other hand we see that 
for an infinitely sized atomic ensemble all m-modes will 
contribute. This is essentially the limit considered in the 
one-dimensional theory in Ref. [Ill ] . That theory applies 
to an infinitely wide sample such that all modes experi- 
ence the same dynamics. For a sample of finite width we 
see that the oscillating behavior of the Bessel functions 
J m gives a cut of the width of the beam scaling with ap- 
proximately r c /<7j_ ~ 1/ \/2J- or r c ~ L/2. This cut r c 
will, due to the behavior of the Bessel function J m , in- 
crease as to increases. We thus see that even though the 
width of the beam is mainly determined by the length 
of the atomic ensemble, the width of the atomic ensem- 
ble plays an important role as a wider ensemble supports 
higher order modes that are inherently wider, thus in ef- 
fect a wider atomic ensemble will generate a wider beam. 

From the expansion Eq. (|VI.4[) and the small argu- 
ment limit of the modified Bessel function Eq. (|VI.9|) 
along with Eq. (|VI.8|) we see that as time increases the 
contributions to the intensity from modes to ^ will 
not grow as rapidly as m = 0. In Fig. [5] we show a 
plot of the radiated power in three SRS modes at time 
t — 0, where we use a Fresnel number T — 4 and optical 
depth d = 160. In Fig. [6] we use an atomic ensemble 
with Fresnel number T = 8 and optical depth d = 160. 
The plots demonstrates how the relative importance be- 
tween different modes are changed as the Fresnel number 
is changed. From the two plots in Figs. [5] and [6] we see 
that the larger the Fresnel number, the more modes with 
higher azimuthal quantum number m can we fit into the 
system. In Fig. [5] we see that the principal mode m = 
is dominating the higher order modes. When the Fresnel 
number is doubled in Fig. [6] the principal mode m = 
is still dominating, but less than in Fig. [5j To conclude 
that a higher Fresnel number, allows higher order az- 
imuthal quantum numbers to to contribute, we have to 
look at the total number of photons for each to. This is 
the topic of Sec. IVI CI and from the results derived there 
we indeed find that we can have relatively more photons 
for higher order to as the Fresnel number is increased. 
E.g for T = 4 the photon power in each of the to = ±1 
mode relative to the to = mode is about 62%, and it 
is 35% for to ~ ±2, whereas for T = 8 this number is 
increased to 72% for to = ±1 and 49% for to = ±2. 

Next we consider how the time evolution changes the 
shape of the mode. From the earlier discussion of Eq. 
(|VI.8j) we expect that the relative photon intensity car- 
ried by modes with to different from the principal mode 




r/a, 



Figure 5: (Color online) Plot of the radiated power for differ- 
ent azimuthal quantum numbers m — 0, ±1, ±2 as a function 
of the detection coordinate r / a± . The plot is taken at the ini- 
tial time, Ft = with an optical depth d = 160. Comparison 
with Fig. [5] demonstrate how the relative distribution of radi- 
ation with different azimuthal quantum number m is changed 
as the Fresnel number T is varied. Here we use 3- = 4 and in 
Fig. we use T — 8. The solid line correspond to m = 0, 
the dashed line correspond to m = ±1 and finally the dotted 
line correspond to m = ±2. 
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Figure 6: (Color online) Same as Fig. [5] but with Fresnel 
number T — 8 



to = will decrease compared to the principal mode as 
time is increased. In Figs. [7] and [8] we plot the radial 
distribution of the photon power at time Tt = 0.25. We 
see that the radial shape of the modes have not changed 
compared with the plots at t = 0, [ Figs. [5] and [6]]. The 
relative maximum photon power for modes with to ^ 
has however decreased compared with the principal mode 
to = 0. Again we can look at the total photon power in 
each mode, and find that for the case of Fresnel number 
T — 4, each of the modes m = ±1 now only contains 
23% of the intensity carried by the to = mode, and the 
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Figure 7: (Color online) Plot of the radiated power for differ- 
ent azimuthal quantum numbers m = (solid line), m = ±1 
(dashed line), and m — ±2 (dotted line) as a function of the 
detection coordinate r/a±. Here the plot is made for a time 
of Ft — 0.25 and an optical depth d = 160. Comparison with 
Fig. [8] demonstrate how the relative distribution of radiation 
in modes with different m is changed as the Fresnel number 
T is varied. Here we use T = 4 and in Fig. [5] we use T = 8. 
When the plots in Figs. [7J and [8] are compared with the plots 
for t — in Figs. [S] and HJ] we indeed see that as time in- 
creases, the evolution of the principal mode, m = is faster 
than that of the higher order modes. 




Figure 8: (Color online) Same as Fig. [7] but with Fresnel 
number T — 8 



to = ±2 mode only 4.8%. A similar behavior is found for 
the T = 8 case, though less pronounced, i.e. now each 
of the modes m — ±1 carries 38% of the photon power 
compared with the m = mode, and for the m = ±2 
modes it is 12%. As expected the modes with m ^ 
become relatively less important for long times. 



C. Total coherent radiation. 



Finally we will examine the total intensity of SRS. We 
shall in this section not only show the effect of the analyt- 
ical calculations made so far, but also compare the result 
with a purely numerical treatment of the equations given 
in Eq. (|II.21|) . The total intensity is normalized such 
that it gives the number of photons per second coming 
through the detector-plane 
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To find the total intensity we use the result in Eq. (|VI.4j) 
and perform the radial integral. To do this we use the 
relation 



[°° rdrJ m (xr)J m (x'r) = 6{X X '\ (VI.11) 
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derived in Appendix [C] The total radiation is then found 
to be 
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Figure 9: (Color online) Plot of the total radiated power P 
measured in number of photons, N p per dr. We in addition 
scale out the natural decay e" rt . With this scaling we get 
a universal curve applying to all ensembles with the same 
Fresnel number. The time axis is scaled in units of dT. We 
use a Fresnel number of T — 4 and show results for three 
different m- modes, m = (solid line), m = ±1 (dashed line), 
and m = ±2 (dotted line). We see that the principal mode 
m — has a slightly faster growth than higher order modes. 

In Fig. [5] we show a plot of the total radiated power, Eq. 
(|VI.12[) for the parameter T — 4. The scaling is chosen 
such that the curve will be identical for all samples with 
the same Fresnel number T . It is interesting to note that 
indeed the intensity in modes with m ^ evolves slower 
in time than for the m = mode. This can be seen by 
looking at the slope of the curves as they are plotted on 
a logarithmic scale. 

In Fig. [TO] we analyze how the total radiation depends 
on the Fresnel number. For large times, the dependence 
is approximately linear in the Fresnel number. This may 
also be concluded directly from Eq. (|VI.12[) . 

We now compare the result for the total radiated power 
with the effective one-dimensional calculation derived in 
Ref. [ll| . The general assumption in the one-dimensional 
calculation is that the atomic ensemble is infinitely wide. 
This assumption makes the problem easy to solve in 
Fourier space. When the transverse momentum in the 
propagator for the light is neglected, the result for the to- 
tal radiated power is that all modes corresponding to dif- 
ferent transverse momentum gives equal contribution to 
the total radiated power. Thus the total radiated power 
measured in units of number of photons per time gives 




ki 

(VI.13) 

Since we have neglected all information on the transverse 
shape there is a priori no upper limit on the transverse 
momentum. Thus taking all modes corresponding to all 
transverse momentum into account gives an infinite con- 




dTt 

Figure 10: (Color online) Plot of the total radiated power 
calculated for varying Fresnel numbers. The lines are calcu- 
lated using the expression Eq. (|VI.12|) for the principal mode 
m = 0. Apart from a complicated behavior at short times 
we see that the total radiation is linearly proportional to the 
Fresnel number. This can also be seen from Eq. (|VI,12[) . The 
solid line correspond to T = 1, the dashed line to T — 2, the 
dotted line to T — 3 and the dash-dotted curve correspond 
to T = 4. 



tribution. A derivation of such a mode des crip tion can 
be found in Ref. [9(. It is concluded in Ref. [12J that for 
a Fresnel number near unity the radiation is dominated 
by a single transverse mode, and thus the total radiation 
is finite, and given approximately by a single term in the 
sum (IVI.13|) . 

We can also make a simplification of our result (|VI.12|) 
by neglecting all kinds of finite size effects in the eigen- 
value matrix, M|,™? n ,. From the derivation of Eq. 
(|VI.12p . one sees that this amounts to fixing {q, q' , I, I'} = 
and setting k = and k' = in the modified Bessel 
function as well as the exponential function. Finally the 
approximation gives an additional factor of 1 + k + k' . 
This is an oversimplification, but allows a comparison 
with the results by Raymer and Mostowski in Ref. [ill ]. 
The total radiated power is then given by 




(VI. 14) 



For f k 1 this expression is identical to a single term 
in the sum in Eq. (|VI.13jl . We now assume the Fres- 
nel number T ~ 1, and apply the approximation (|VI.9|) . 
which is only valid for small Fresnel numbers. In this 
way we find 




(VI. 15) 
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and the integral results in the total radiated power 
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We are thus led to conclude that for a Fresnel number 
near unity, the simple Raymer Mostowski result corre- 
spond to neglecting all spatial corrections to the dynamic 
of the atoms and also neglecting spatial corrections to the 
propagation of light out of the atomic ensemble. 

We can improve the approximation, by looking at the 
general result in Eq. (|VI.12[) and keeping only zeroth 
order terms in the index q, q' , I and In this way we get 
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for T <C i + we can reduce Eq. (|VI.17|) even further 
and arrive at the result 
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In this limit T <C \ + ^p- the only contribution to the 
total radiated power comes from the m = mode. 

In Fig. [IT] we analyze how the different corrections to 
the Raymer Mostowski calculation effects the total ra- 
diated power. We fix the Fresnel number at T — 1, as 
this is the limit where the Raymer Mostowski result is 
assumed to be valid. The curve Po(t) is the simple result 
Eq. (|VI.16[) . In curve P\ (t) we use the lowest order finite 
size correction, that is Eq. (|VI.17[) . Finally in curve P(t) 
we use the general result from Eq. (|VI.12[) . which is eval- 
uated nummerically with the approximated Bessel func- 
tion (|VI.9|) . The approximation is used so that we gan 
get an estimation of the effect of all azimuthal quantum 
numbers, thus the nummcrical methods require modest 
Fresnel numbers, such that the main contribution to the 
total radiated power comes from the mode coresponding 
to m = 0. We see that the simple Raymer Mostowski 
type result, Eq. (|VI.16|) over-estimates the total radi- 
ated power compared to the general result. We also see 
that the zeroth order result P\(t) is a much better ap- 
proximation in the regime dTt/8 3> T . 

Finally we compare the result of Eq. (|VI.12[) with a 
purely numerical calculation based on the point particle 
equations (|II.21[) and (|II.22[) . To make such a comparison 
we need to connect the evolution of the atomic operators 
bj (t) with the total intensity of the radiated field. Based 
on energy conservation, the evolution of the number of 
atoms in the ground state, is given by the number of 
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Figure 11: (Color online) Plot of the total radiated power 
P(t) scaled so that it only depend on Fresnel number. Here 
we use T = 1. To demonstrate the effects of a finite sized 
atomic ensemble, we show three different curves. The solid 
line is P(t), the general result from Eq. (|VI.12[) . The dashed 
line is Pi(t) (|VI.17|) . where we use the zeroth order expansion 
of the general result (|VI.12|) and assume a large value of dTt. 
Finally the dotted line is Po(t) ()VI.16[) . where we completely 
neglect all geometric effects on the matrix M{!!^} n , . 



photons exiting a boundary sphere enclosing the atomic 
ensemble. We derive this conservation law in Appendix 
ID] where we show that 
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where Mjj> is given by Mjj> + FSjji , and Mjji is given in 
Eq. (jll.lfp . When comparing the result of Eq. (|VI.12p 
to the atomic evolution we have to remember that we 
are only measuring half of the photons, since we only 
consider the emission at one end of the ensemble. Using 
that the evolution of the atomic operators are given by 
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we find that the atomic operators evolve in time accord- 
ing to 



(VI.21) 



where we define M as the matrix with elements given by 
Mjj> . After taking quantum average of the result in Eq. 
(|VI.19|) we find that 



k s he 



dQ (D ( + ) {r,t)D^>(r,t)) 



(+)/ 



trace e 



M't 



{t + M)e Mt \ +C.c 



(VI.22) 



17 



We then find the total intensity from the point particle 
model 

PnQ) = \ {trace [e £rtT {t + M)e &t ) + C.c.}, (VI.23) 

where we normalize with a factor 1/2 since we want to 
compare the result with the result in Eq. (|VI.12j) . 

The advantages of making these calculations, or indeed 
solving the problem of SRS on a computer are clear. One 
avoids the problems of shifting from the point particle 
model to a continuous model. Thereby one also automat- 
ically include dipole dipole interaction effects connected 
to the point particle nature of the system which we have 
ignored here. Also the computer easily describes the total 
radiated field and not only the strongest super-radiating 
mode as we have analyzed here. On the other hand the 
direct method is numerically heavy for a large number 
of atoms, and we are limited to N ~ 6000 atoms. To 
understand the behavior at larger number of atoms it is 
therefore important to have an analytical theory along 
the lines considered here. 

To make the numerical simulation we have randomly 
distributed between 3000 and 6000 atoms with a distribu- 
tion function given by Eq. (|III.2|1 . After that the matrix 
Mjji is calculated and processed in order to find the to- 
tal number of Stokes photons (|VI.23[) . We can then by 
making a series of such realizations of the position of the 
atoms get some statistics on the inherent noise on the 
point particle model. In Fig. [12] we show the result of 
a numerical calculation using parameters T = 4 and an 
optical depth of d — 90. When we increase the num- 
ber of atoms, we decrease the particle density in order to 
keep a fixed Fresnel number and a fixed optical depth. 
We see from Fig. [T2] that there is some dependence on 
particle density, an effect of the fact that the system is 
a point particle system and not a continuum, hence we 
do not expect the analytical theory developed so far to 
explain this effect. However as the density desceases the 
total radiated power converges. Finally in Fig. [T3] we 
compare the total radiated power in the analytical cal- 
culation P(t), (|VI.12[) . with the nummerical calculation 
P N {t), (|VI.23|) . That the two methods gives very dif- 
ferent results for small times is quite clear since initially 
the radiation is dominated by the spontaneous emission, 
which is not included in the analytical calculation. At 
increasing times, which is the regime where the analyti- 
cal calculation is supposed to be valid, the two methods 
gives quite similar results, and we therefore believe that 
the analytical calculation gives an accurate description. 

We finally note that for the time-scale used in Fig. 
1131 the approximation of neglecting depletion is not 
completely justified, as the number of emitted photons 
exceeds the number of atoms already before the two 
curves meet. We can examine the break-down of the no- 
depletion assumption, by finding the time t c , at which the 
number of photons emitted in the superradiating mode 
Np(t) exceeds the number of atoms in the ensemble Na, 




Figure 12: (Color online) Nummerical results for the total 
radiated power Pjv(£) per decay rate T in the point particle 
model (IVI.23|l . To exploit the nummerical model we fix the 
Fresnel number T = 4 and the optical depth d = 90, but vary 
the number of atoms involved. The solid line correspond to 
TV = 6000 atoms, the dashed line to N = 5000, the dotted 
line to TV = 4000, and finally the dash-dotted line correspond 
to TV = 3000 atoms. As the plot shows there is a dependence 
on the atomic density due to point particle effects that is not 
included in the analytical theory, but as the atomic density 
decreases (with increasing TV) the curves seem to converge. 
The errorbars indicate the noise inherent in the point particle 
model due to the random positions of the atoms. 




Figure 13: (Color online) Comparison of the analytical calcu- 
lation of the radiated power P(t), (|VI.12|) . (solid line), with 
the nummerical result P/v(i), (|VI.23[) . (dashed line). The 
Fresnel number is T = 4 and the optical depth is d — 90. 



i.e. Na — Np(t c ), where 



N P (t) = / P{t')df 
Jo 
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To get an analytical result we will use the approximation 
P(t) « Pi(t), with Pi{t) given in Eq. (|VI.17p . After the 
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Figure 14: Total number of photons emitted in the SRS mode 
N p (t) divided by Fresnel number T as a function of the scaled 
time dFt. 



integration in Eq. (|VI.24[) we find 
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where we have used that <1,Na 3> 1. In Fig. [14] we plot 
the function N p (T). From the requirement Na = N p 
we find the time t c for the result shown in Fig. [13] to 
be Tt c = 0.54, where the radiation is still dominated by 
spontaneous emission. If, however we increase the optical 
depth we decrease the time at which the analytical curve 
P(t), (|VI.12|) . and the numerical curve P N (t), (|VI.23|> . 
agree. For a higher optical depth d there will thus be 
a region where the effects considered here are dominant 
within the applicability of our theory. While the limited 
atom number used here is thus not physically relevant, 
the simulation can still be used as a confirmation of the 
approximations used in our analytical calculation since 
both curves are derived using the same approximation 
of neglecting depletion of the atoms. For the ongoing 
SRS experiments using Bose-Einstein condensed atoms 
e.g. Ref. [i[ the number of atoms used in the process 
is factors of thousands larger than what we are able to 
numerically simulate here, and the approximation used 
here is much less severe. 



VII. CONCLUSION 

In this paper we have developed a three-dimensional 
theory for spontaneous Raman scattering (SRS). The 
theory applies to an ensemble of non-moving atoms and is 
derived by describing the atoms as a continuous medium. 
In the theory we neglect the depletion of the initial 
atomic state and the theory is therefore mainly applicable 



to the onset and build up of SRS. We believe, however, 
that the theory still captures the most important effect 
of the three-dimensional structure of the problem, since 
after the onset of SRS the radiation is dominated by the 
modes determined by our theory. 

The theory is based on a generalization of the one- 
dimensional theory in Ref. [111 ]. In the limit where the 
Fresnel T is very large we find that the one-dimensional 
description of Ref. [llj applies to all transverse modes 
in agreement with the derivation in Ref. [9(. Without a 
detailed investigation of the three-dimensional structure 
there is, however, no restriction on the transverse mo- 
mentum of the light and a naive application of the theory 
therefore predicts an infinite radiated intensity. In our 
three-dimensional theory the build up of SRS only hap- 
pens for small transverse momentum of the light, limited 
by the Fresnel number T of the ensemble. This automat- 
ically limits the emitted radiation such that the theory 
gives finite predictions. 

In the theory we assume that the dimensions of the 
ensemble is much larger than the wavelength of the ra- 
diation, and we show that in this limit the only two pa- 
rameters describing SRS are the optical depth d and the 
Fresnel number F of the ensemble. We find that in the 
limit T > 1 the time scale of SRS is almost exclusively 
given by the optical depth d with only a weak depen- 
dence on the Fresnel number T . On the other hand, the 
total radiated power for SRS depends strongly on the 
Fresnel number T . The total power radiated into modes 
with a given azimuthal quantum number m is linearly 
proportional to T . We find that the largest contribu- 
tion to the radiation always comes from the azimuthal 
quantum number m = and that this contribution also 
have the fastest growth. With increasing Fresnel number 
T the contribution from other azimuthal quantum num- 
bers m 7^ may, however, become comparable to the 
m = contribution. To investigate the validity of our 
analytical findings we have compared our analytical re- 
sults to a direct numerical solution for a limited number 
of atoms. The two approaches are found to be in good 
agreement. 

An interesting question which we have not addressed 
in detail here comes from the fact that an ensemble of 
atom is not given by a continuous density, but consists 
of a collection of discrete point particles. The effect of 
this is in principle included in our direct numerical inves- 
tigations, and may be the reason for the dependence on 
the atom number in Fig. [T2] Here the simulations with 
the highest density deviate from the results with a lower 
density. It would be interesting to investigate such effects 
using for instance the methods developed in Ref. [lj|. 
Furthermore, the question of the collective emission of 
radiation from atomic ensembles is also very interesting 
from the point of view of quantum information. Several 
important quantum protocols such as quantum repeaters 
26 1 , quantum memory [2?| , and quantum teleportation 
28| are currently being investigated in atomic ensembles. 



For a full evaluation of the potential of these approaches 
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it will by important to have a full understanding of the 
effect of the realistic three dimensional structure of the 
ensembles. The methods developed in this article may 
serve as useful starting point for such investigations. 
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Appendix A: DERIVING THE FIRST ORDER 
CORRECTION TO THE MATRIX M^„, 



Equation (|A.1|) may be rewritten as 
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From Ref. [24| we find the integral to give 
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By introducing the dummy variable a = 2a\ in the 
Gaussian function, the series expansion of the x-integral 
in Eq. (llV.llj) may be written as 
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We see that in terms of an expansion in the variable l/oj_ 
each differentiation will give a factor of l/oj_. We shall 
therefore only consider a sum up to the first order in the 
differential. To zeroth order the x-integral simply gives 



Using the above expansion along with the relation 
Im{x) = i~ m Jm{i x ) together with the result [24| 
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To first order we find the x-integral to give 
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To understand the above expression let us assume a sufficiently large a± so that the modified Bessel function I m ±i 
can be approximated with I m . In this way we get 
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The above approximation gets worse for increasing values width of the sample, higher order modes in m has less 
of to, however we argue in Sec. IVI B[ that for a finite influence. Finally the exponential function along with 
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the modified Bessel function express a conservation of 
transverse momentum given by the variables 7„ since for 
increasing values of the transverse momentum, Eq. (|A.7j) 
can be approximated with 
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We shall then make the approximation 
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thus the expression in Eq. (|A.8|) can to second order in 
the difference 7„ — 7„< be written as 
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This result is the large size limit, and we therefore con- 
clude that to give this limit as a± — ► 00 the term in Eq. 
(|A.6|) must be approximated with 
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From this we conclude the result given in Eq. pV.13p . 



Appendix B: COMMUTATION RELATION FOR 

A™„, AND A 1 ™ , 

Here we show that the two matrices A™ n , and A 1 ^, 
commute. Since both matrices are symmetric, it is 
enough to show that the product J2 p AJ^A 1 ^, is sym- 
metric. Again we make the continuation Y\ — > f 

G _ t-~<p a c J 7r 

for a c — > 00. In this way we get 
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After making the 7 P -integral we end up with 

4^4 -4(72+7^)/ ~ ,) 

2_^A™ A 1 ^, = — ^-j-p ^ — — — — ^ — ■ (B.2) 
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Since the matrix Eq. (|B.2|) is symmetric we conclude 
that the matrices A™ / and A 1 ™/ commute. 



Appendix C: DERIVATION OF EQ. (|VI.11[) 



Here we will show Eq. (|VI. 1 1|) . Our starting point is 
the orthogonality relation given by 

/ rdr Jm(n n r)J m {ln>r) = — ""' c (C.l) 



where the variable j n = mn and X rnn is the n'th zero of 
the m'th order Bessel function J m . We will assume that 
X mn is large, which does not require 7 n to be so, since 
we can choose the cut-off a c to be anything. In this way 
we can write Eq. (jC.ip as 



rdr J m h n r)J m (ln'r) = - 

^In 



(C.2) 



We will then take the sum over n on both sides and use 
the standard continuation Y7, — > f ^ 22L so that 



d^nln I rdrJ m ("f n r)J m (~f n ,r) = 1. (C.3) 



Since 7™ is now a continuous variable, we conclude that 
the measure of the distribution 

f(x,x') = x I rdrJ m (xr)J m (x'r), (C.4) 



where x, x' is some real and positive number is unity. 
The next step is to show that for x 7^ x' the function 
f(x,x') vanish. This follows when choosing a zero point 
X mn and a cut-off a c such that say x = 7„. This does 
not necessarily mean that x' has a similar representation 
with the chosen cut-off. On the other hand this is not 
necessary as one may show, see e.g. [2l|, that 

(Cf 2 n -x ) / rdrJ m (j n r)J m (x'r) = 0. (C.5) 
Jo 

from here we conclude that when 7„ and x 1 are different 
the function f(-~f n ,x') vanish. This concludes the deriva- 
tion of Eq. ifVLTT]) . 



Appendix D: THE SUM RULE 



Here we derive the sum rule Eq. (|VI.19|I used in Sec. 
IVI CI The starting point is the total radiated intensity 
of Stokes-photons 

j> |d _ x (V x A+) - (V x A") x D+j, (D.l) 

where S is a sphere surrounding the atoms. Using the 
Divergence theorem as well as the Maxwell equations, the 
total radiated intensity can be written as 
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To lowest order in 1/uj s , Eq. (|D.2|1 reduce to 
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where we have used Eqs. (|II.11|) . (|II.17[1 and (|II.18|) . 

When measuring the intensity infinitely far away from 
the atomic ensemble, the expression in Eq. (jD.lj) reduce 
to the electric field squared times 2/ioc, thus the normal- 
ized sum-rule reads 

^-/dQD-.D + = ^rS, W St W 

+ ^2{b j (t)M jj ,b],{t) + H.c.}. (D.5) 
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